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Large-eddy simulation of 
flow past a circular cylinder 

By R. Mittal 

1. Motivation and objectives 

Some of the most challenging applications of large-eddy simulation are those in 
complex geometries where spectral methods are of limited use. For such applications 
more conventional methods such as finite difference or finite element have to be used. 
However, it has become clear in recent years that dissipative numerical schemes 
which are routinely used in viscous flow simulations are not good candidates for 
use in LES of turbulent flows. Except in cases where the flow is extremely well 
resolved, it has been found that upwind schemes tend to damp out a significant 
portion of the small scales that can be resolved on the grid. Furthermore, it has 
been found that even specially designed higher-order upwind schemes that have 
been used successfully in the direct numerical simulation of turbulent flows produce 
too much dissipation when used in conjunction with laxge-eddy simulation. 

A case in point is the LES of flow past a circular cylinder performed by Beau- 
dan &; Moin (1994) at a Reynolds number of 3900. One of the objectives of this 
investigation was to study the suitability of higher order upwind-biased schemes 
for LES of complex flows and to validate the methodology against experimental 
results of Ong & Wallace (1994) and Lourenco & Shih (1993). In particular, 5 th - 
and 7 th -order schemes were used for these simulations. The 5 tA -order scheme has 
been successfully used for DNS of transition and turbulence in flow over a flat plate 
by Rai & Moin (1993) and it was thought that these schemes would be useful in 
LES of flows in complex geometries. However, the conclusion of the study by Beau- 
dan Sz Moin (1994) was that except in regions where the mesh was fine enough to 
resolve a significant portion of the small scales, numerical dissipation overwhelmed 
the contributions from the subgrid-scale eddy- viscosity model. 

In contrast to upwind-biased schemes which control aliasing through numerical 
dissipation, aliasing is controlled in central schemes by an energy conservation prin- 
ciple. Such schemes do not exhibit numerical dissipation and, therefore, there is 
no spurious damping of the smaller scales. This feature makes the schemes attrac- 
tive for use in LES of complex flows. The downside of using such schemes is the 
dominance of dispersive error, which makes these schemes extremely sensitive to 
aspects such as the grid stretching factors (Cain & Bush, 1994) and outflow bound- 
ary conditions (Gresho & Lee 1981). Thus, even though the central schemes might 
have a clear advantage over upwind biased schemes in simple geometries, in more 
complex geometries where complicated grids are used, the superiority of central 
schemes needs to be established and this is the motivation of the current study. 

The objective of the current study is to perform a LES of incompressible flow past 
a circular cylinder at a Reynolds number of 3900 using a solver which employs an 
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energy- conservative second-order central difference scheme for spatial discretization 
and compare the results obtained with those of Beaudan &: Moin (1994) and with 
the experiments in order to assess the performance of the central scheme for this 
relatively complex geometry. 

Salient features of the simulation of Beaudan & Moin (1994): 

Beaudan & Moin (1994), henceforth refered to as BM, simulated the flow past a 
circular cylinder on an 0-mesh using a compressible flow solver. One point upwind- 
biased 5 th - and 7 -order schemes were used for the spatial discretization of the 
convective terms. The simulations were carried out on a 144x136x48 ( rxOxz ) grid 
and good resolution was provided near the cylinder surface and in the near wake 
region (x/D < 2.0). Beyond this region, the grid was stretched geometrically in 
the streamwise direction such that the streamwise grid spacing at x/D = 10.0 was 
about 0.13H. The mesh near the outflow boundary was made extremely coarse in 
order to damp out disturbances and a convective outflow boundary condition was 
used. Grid stretching ratios in excess of 10% were used to obtain the desired grid 
spacing in the wake. 

Simulations were carried out with no subgrid-scale model, with a fixed coefficient 
Smagorinsky model and with the spanwise averaged version of the dynamic model 
(Ghosal et al., 1995, Moin et al. 1991). It was observed that mean wall statistics 
such as drag, pressure coefficients, wall shear stress and separation angles were not 
significantly different in the three simulations and all showed reasonable agreement 
with experimental data. In the vortex formation region ( x/D < 4.0), it was found 
that the dynamic model predicted mean velocities and Reynolds stresses which were 
in better agreement with the experimental results than the other two simulations. 
Beyond this region the difference between the three computed solutions diminished 
such that the solutions were virtually indistinguishable beyond x/D > 7.0. It was 
found that in this region where the mesh was relatively coarse, numerical dissipation 
overwhelmed the contribution of the SGS model. The simulation with the 7 th -order 
scheme showed evidence of increased energy in the high wavenumbers, but here 
too it was found that a substantial portion of the resolvable wavenumber range 
was damped due to numerical dissipation. It was concluded that these high order 
upwind-biased schemes were unsuitable for use in LES. 

2. Accomplishments 

2.1 Numerical method 

The solver used in the current work is based on the solver developed by Choi et 
al. (1992) and employs a second-order central-difference method written in general- 
ized coordinates in a spanwise periodic domain. Velocity components and pressure 
axe fully staggered in order to strictly conserve mass in the generalized coordinates. 
It should be pointed out that strict conservation of momentum and energy is not 
guaranteed on a non-equispaced mesh. The solution is advanced in time using a frac- 
tional step scheme wherein a third-order Runge-Kutta scheme and a Crank-Nicolson 
scheme is used for the nonlinear convection terms and viscous terms respectively. 
A multigrid solver is used in conjunction with a Gauss-Siedel line-zebra scheme 
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FIGURE 1. C-mesh used for Run-II. 


for solving the pressure Poisson equation. The solver employs a spanwise- averaged 
version of the dynamic model where the total viscosity is constrained to be greater 
than zero (Ghosal et al ., 1995). The spanwise length of the cylinder is chosen to be 
7 tD which is the same as BM. 

A C-mesh is used for the present simulations (Fig. 1). This type of mesh is ideally 
suited for simulating wake flows since better streamwise resolution can be selectively 
provided in the wake region. The use of a C-mesh also simplifies the application of 
outflow boundary conditions. Furthermore, another advantage of using a C-mesh 
is that as the flow separates from the cylinder, it remains roughly aligned with one 
family of grid lines, and thus good control over the streamwise stretching ratio can 
be maintained in this region. It has been found that in LES, where the resolution 
is at best marginal, central schemes can tolerate only a small streamwise stretching 
factor (< 3%). Higher stretching factors can leads to the amplification of grid-to- 
grid oscillations (2 — A waves). If an O-type mesh were to be used for the present 
simulations, the flow in the region of the separated shear layer would experience 
large strething ratios as it would go from being aligned with one family of grid lines 
to being aligned with the other, and solution in this region would be contaminated 
by 2 — A waves. Thus, the use of a C-mesh is necessary for obtaining a good solution 
with the current solver. This brings in the important point that the grid has to be 
designed keeping in mind the underlying spatial discretization. 
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2.2 Simulation results and discussion 

The first simulation (Run I) was carried out on a 329x100x48 mesh with 80 
points on the cylinder surface, 125 streamwise points along the wake centerline, 
100 points in the wall normal direction, and 48 points along the spanwise direction. 
Since this was the first simulation, a relatively coarse mesh was chosen with the 
objective that results from this simulation would provide an estimate of the reso- 
lution requirements. The results from this simulation are summarized in Table 1. 
Also tabulated for direct comparison are the corresponding results from the 2-D 
simulation and 3-D LES of BM and experimental results from various studies. This 
simulation predicted a higher mean drag, rms lift, and base suction pressure coef- 
ficient than the corresponding LES of BM and experiments. Furthermore, it was 
observed that the computed in-plane Reynolds stresses (u' 2 , v' 2 and u'v') in the near 
wake were significantly higher than the corresponding LES of BM and experiments 
of Lourenco & Shih (1993). On the other hand, spanwise Reynolds normal stress 
(uP 2 ) in the near wake was under-predicted. All indications were that the flow was 
not developing enough three-dimensionality. 

To get a realistic evolution of the three-dimensionality in the near wake, one re- 
quires adequate resolution of the underlying two-dimensional flow in addition to 
good spanwise resolution of the three-dimensional structures. It was clear that the 
azimuthal resolution of the attached boundary layer and separation region was much 
less than in the LES of BM. This could possibly lead to an incorrect location of 
the separation point and subsequent evolution of the separated shear layer. There- 
fore, it was decided to continue the simulation on a mesh with increased azimuthal 
resolution on the cylinder surface. 

The second simulation (Run-II) was carried out on a 399x100x48 mesh where 
the number of points on the cylinder surface was increased from 80 to 150. In order 
to maintain a smooth streamwise distribution of grid points at the concave corner 
in the base of the cylinder, streamwise resolution had be improved marginally (by 
about 10%) in the near wake. The grid was kept roughly the same in all other 
regions. Some of the results of this simulation are summarized in Table 1. It was 
observed that overall there was no substantial improvement in the results. The 
mean drag coefficient, rms lift coefficient, and base pressure coefficient all show a 
small change towards the correct values but the results are still significantly different 
from BM and experiments. 

In Fig. 2 is shown the variation of lift and drag coefficient with time after the 
flow has reached a statistically stationary state. All the data presented for this 
simulation has been averaged over the time period shown in this figure. Figure 3 
shows the distribution of the surface pressure coefficient obtained from the present 
simulations. Results of BM have also been plotted for comparison. It is clear that 
the current simulations predict a significantly higher suction pressure in the wake 
region and that increasing the azimuthal resolution on the cylinder surface has only 
a marginal effect on the surface pressure distribution. 

Figures 4a and 4b show the streamwise and cross-stream mean velocity profiles 
in the near wake (x/D = 1.54) obtained from the current simulations. Figure 4a 
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Table 1. Wall Statistics 



Run — 

I Run — II 

BeaudanSzMoin 

(3-D) 

BtaudanfzMoin 
( 2-D) 

Lourenco 

hShih 

Cp, 

-1.28 

-1.15 

-0.95 

-2.16 

— 0.9±0.05 

Cd 

1.2 

1.1 

1.0 

1.74 

0.98±0.05 

0s 

00 

CO 

© 

© 

00 

00 

85.8° 

108.1° 

85°±2° 



FIGURE 2. Variation of lift and drag coefficient with time obtained from Run-II. 

Cd\ Cl- 

shows that the current simulations underpredict the momentum deficit in the near 
wake, and consequently the wake bubble length is also underpredicted (see Table 
1). In contrast to the streamwise velocity, the mean cross-stream velocity (Fig. 4b) 
matches well with the results of BM. Furthermore, it is observed that the experi- 
mental data does not match with any of the simulation results. This is consistent 
with the fact that Beaudan & Moin (1994) indicated that large errors might be 
present in the experimental measurements (Lourenco & Shih, 1993) of cross-stream 
velocity in the near wake. 
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FIGURE 4. Mean velocity profiles at x/D = 1.54. (a) Streamwise velocity, (b) 

Cross-stream velocity. : Run-II; : Run-I; • : Beaudan & Moin; * : 

Lourenco & Shih. 

In Fig. 5 are shown Reynolds stress profiles at this streamwise location. It can 
be observed from Fig. 5a that the current simulations over-predict the streamwise 
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Figure 5. Reynolds stress profiles at x/D = 1.54. (a) u n (b) v' 2 (c) u f v t2 (d) w* 2 
: Run-II; : Run-I; • : Beaudan & Moin; * : Lourenco & Shih. 

Reynolds stress (u' 2 ). However, overall, the stress profile is better predicted in Run- 
II. The noticeable asymmetry of the profile about the wake ceneterline obtained 
from Run I also suggests that more than six shedding cycles might be needed for 
averaging the statistics. Figure 5b shows the corresponding profiles of cross-stream 
normal Reynolds stress (v' 2 ), and here large differences between the results of the 
current simulations and the results of BM can be seen. Run-I and Run-II over- 
predict the peak stress by about 100% and 80% respectively. A similar trend is 
observed in Fig. 5c, in which profiles of (u , v , ) axe plotted. 

Figure 5d shows profiles of the spanwise Reynolds normal stress, and we observe 
that the current simulation under-predicts this stress component. There is however, 
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FIGURE 6. One-dimensional spectra at x/D = 0.7, y = 0. :Run-II; • : 

Beaudan & Moin. 


a noticeable effect of increasing the azimuthal surface resolution on this stress com- 
ponent. First the peak stress obtained from Run-II is about 17% higher than that 
obtained from Run-I. A more noticeable effect is that in Run-II a ‘fuller’ profile is 
obtained in the region y/D>0.25 and this is in much better agreement with BM. 
This is most likely due to improved streamwise resolution in the near wake which 
leads to increase in the growth of three-dimensional instability in this region. 

Figure 5 shows clear evidence that in the current simulations, the flow is not 
developing enough three-dimensionality in the near wake. As a result of this, in- 
plane stresses are over-predicted and spanwise stresses axe under-predicted. It ha s 
been shown that the in-plane Reynolds stresses play a significant role in determining 
the base suction pressure (Mittal &; Balachandar, 1995). Thus the higher in-plane 
Reynolds stresses lead to a higher base suction pressure and drag in the current 
simulations. 

Figure 6 shows the spanwise one-dimensional spectra of the streamwise velocity in 
the near wake obtained from Run-II and the simulation of BM. Direct comparison 
of the spectra can be made since both simulations employ the same resolution in 
the spanwise direction. It can be observed that the two spectra match well only 
for the low wavenumbers (approximately 20% of the wavenumber range). Beyond 
this range, the spectra obtained by BM exhibits significant damping and the energy 
shows a decay of about seven orders of magnitude. In contrast, the spectra obtained 
from the current simulation is relatively flat with about one order of magnitude 
decay in the high wavenumber range. It should be pointed out that comparison of 
spanwise spectra at other wake locations shows a similar trend. Thus, it is clear 
that the higher-order upwind scheme used in the simulation of BM damps out a 
significant portion of the wavenumbers that can be resolved on the grid. 
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3. Summary and future plans 

The results indicate that the evolution of the secondary instability that is respon- 
sible for the generation of three-dimensionality in the near wake is not captured 
well in the current simulations. This discrepancy could result from inadequate res- 
olution of the underlying two-dimensional flow and/or spanwise resolution of the 
three-dimensional structures. The current simulations have the same number of 
grid points in the spanwise direction as Beaudan & Moin (1994). However, given 
that the current simulations use only a 2 nd -order accurate spatial discretization, 
more spanwise grid points might be needed to match the resolution power of the 
5 th -order scheme. Comparison of the modified wavenumber for the schemes sug- 
gests that the 2 nrf -order scheme might need up to twice the number of grid points 
to match the resolution of the 5 th -order scheme. 

It is also clear from the present study that the restriction imposed on the stream- 
wise grid stretching factor when using central schemes represents a severe constraint 
on mesh design for complex geometries. In this respect, the higher-order upwind 
biased schemes are more flexible since they allow the use of higer stretching factors 
and increased resolution can be provided selectively at desired locations. However, it 
is also evident that even these higher order upwind schemes exhibit significant dissi- 
pation, and the scales corresponding to the top half of the wavenumber range, which 
are crucial for determining the subgrid-scale dissipation, are effectively damped out 
due to the numerical dissipation. The second-order central difference scheme, on 
the other hand, preserves the energy in the small scales and allows the subgrid-scale 
dissipation to have a more significant impact on the resolvable flow field. 

Doubling the number of grid points on the cylinder surface improves the results 
only marginally. Therefore, it is unlikely that the disagreement in results is due 
to lack of resolution on the cylinder surface. In-plane resolution in the near wake 
region could also be one cause of the discrepancy. In particular, the restriction on 
the streamwise stretching ratio and the presence of the concave corner at base of 
the cylinder result in the near wake having poorer streamwise resolution than the 
simulation of BM. A systematic spanwise resolution study would require doubling 
the spanwise grid points which would effectively double the computational resources 
required. In contrast, doubling the streamwise resolution in the near wake can be 
accomplished with about a 30% increase in computational resources and is thus the 
more viable next step. 

The near-term objective then is to obtain wall and near wake statistics which 
are independent of the near wake in-plane resolution. 2-D simulations, which are 
relatively cheap, can be used to give a rough estimate of the in-plane resolution 
requirement. Once statistics which are independent of the in-plane resolution in 
the near wake are obtained, these will be compared with the results of BM. If the 
wall and near wake statistics match reasonably well with BM, this will imply that 
the spanwise resolution is adequate and the next step will then be to obtain and 
compare the statistics in the downstream wake region. On the other hand, if the 
wall and near wake statistics do not match with BM, this will be an indication 
that increased spanwise resolution might be required. The code is in the process of 
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being ported to the IBM SP2 parallel computer where the turnaround time will be 

significantly reduced and it will be possible to use larger meshes. 
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